import data from 86 files .nc4¶

  • date: 2019-01-01
    • variables:
      • latitudes
      • longitudes
      • CO2_std
      • CO2_sdev
    • dimensions:
      • time
      • lat
      • lon
    • attributes:
      • year
      • month
      • day
      • time

Normailize data ???¶

In [ ]:
# import libs
import os
import netCDF4 as nc4
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
import itertools
In [ ]:
# Ruta al directorio que contiene los archivos netCDF
directorio_archivos = 'descargas'
# Lista de archivos netCDF en el directorio
archivos_netCDF = [os.path.join(directorio_archivos, nombre_archivo) for nombre_archivo in os.listdir(directorio_archivos) if nombre_archivo.endswith('.nc4')]
# create a dict de dimensiones [91,144]
row = np.arange(0,91)
col = np.arange(0,144)
dict_dim = dict()

def media_armonica(valores):
    inversos = [1 / x for x in valores if x != 0]  # Calcula los inversos de los valores no iguales a 0
    if len(inversos) == 0:
        return 0  # Si no hay valores válidos, devuelve 0
    return len(inversos) / sum(inversos)

for i,j in itertools.product(row,col):
    dict_dim[(i,j)] = {"lat": [], "lon": [], "co2_std": [], "co2_sdev": []}
    # dict_dim[(i,j,"lon")] = []
    # dict_dim[(i,j,"co2_std")] = []
    # dict_dim[(i,j,"co2_sdev")] = []

for i in range(len(archivos_netCDF)):
    # Abre el archivo netCDF
    data = nc4.Dataset(archivos_netCDF[i],"r")

    latitude = data.variables['Latitude'][:]
    longitude = data.variables['Longitude'][:]
    co2_std = data.variables['mole_fraction_of_carbon_dioxide_in_free_troposphere'][:]
    co2_sdev = data.variables['mole_fraction_of_carbon_dioxide_in_free_troposphere_sdev'][:]
    # print(latitude)
    # print(longitude)
    # print(co2_std)
    # print(co2_sdev)
    lon, lat = np.meshgrid(longitude, latitude)

    # Aplana las mallas de coordenadas y la matriz de co2_free_troposphere a vectores unidimensionales
    lon_flat = lon #.flatten()
    lat_flat = lat #.flatten()
    # create data.fream
    df_lat = pd.DataFrame(lat_flat).fillna(0)
    df_lon = pd.DataFrame(lon_flat).fillna(0)
    df_co2_std = pd.DataFrame(co2_std).fillna(0)
    df_co2_sdev = pd.DataFrame(co2_sdev).fillna(0)




    for i,j in itertools.product(row,col):
        aux_lat = df_lat.iloc[i,j]
        dict_dim[(i,j)]["lat"].append(aux_lat)
        dict_dim[(i,j)]["lon"].append(df_lon.iloc[i,j])
        dict_dim[(i,j)]["co2_std"].append(df_co2_std.iloc[i,j])
        dict_dim[(i,j)]["co2_sdev"].append(df_co2_sdev.iloc[i,j])     

# Crear un data frame de los puntos unicos lat, lon para la dimension i,j
        
# print(dict_dim[(0,0)]["lat"])
columns = ["co2_std", "co2_sdev", "lat", "lon"]
df_dim = pd.DataFrame(columns=columns)
glat, glon = [], []
for i,j in itertools.product(row,col):
    # print(dict_dim[(i,j)]["lat"])
    # print(dict_dim[(i,j)]["lon"])
    # print(dict_dim[(i,j)]["co2_std"])
    # print(dict_dim[(i,j)]["co2_sdev"])
    # print("==================================")
    # print("==================================")
    # print("==================================")
    # print("=============================
    x = media_armonica(dict_dim[(i,j)]["co2_sdev"])
    y = media_armonica(dict_dim[(i,j)]["co2_std"])
    #print(pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0])
    glat.append(pd.Series(dict_dim[(i,j)]["lat"]).unique()[0]), glon.append(pd.Series(dict_dim[(i,j)]["lon"]).unique()[0])
    min_lat_lon = min([len(pd.Series(dict_dim[(i,j)]["lat"]).unique()),len(pd.Series(dict_dim[(i,j)]["lon"]).unique())])
    if min_lat_lon == 1:
        df_dim.loc[len(glat)-1] = [x, y, pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]
        # ({"co2_std": x, "co2_sdev": y, "lat": [pd.Series(dict_dim[(i,j)]["lat"]).unique()[0]], "lon": [pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]})
    else:
        df_dim.loc[len(glat)-1] = [x, y, pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]
        print("no es unico".upper().ljust(30,"="))
    # print()    
                

    
        # value.append(df_lat.iloc[i,j])
        # value.append(df_lon.iloc[i,j])
In [ ]:
# normalize df_dim
x = range(6)
np.std(x)
np.mean(x)
Out[ ]:
2.5
In [ ]:
df_dim["co2_sdev"]
Out[ ]:
0        0.000395
1        0.000407
2        0.000401
3        0.000394
4        0.000403
           ...   
13099    0.000000
13100    0.000000
13101    0.000000
13102    0.000000
13103    0.000000
Name: co2_sdev, Length: 13104, dtype: float64
In [ ]:
#create a geodataframe with geopandas and shapely: lat, lon -> glat, glon 
#from matplotlib.patches import Polygon
#geometry = [gdp.Polygon(xy) for xy in zip(glon, glat)]
geometry = [Point(xy) for xy in zip(glon, glat)]
crs = 'EPSG:4326'
gdf = gpd.GeoDataFrame(df_dim, crs=crs, geometry=geometry)
In [ ]:
min(df_dim["co2_std"]) 
float(max(df_dim["co2_std"]))
df_dim.query("co2_std >= 0.000005")
Out[ ]:
co2_std co2_sdev lat lon
269 0.000005 0.000396 88.0 132.5
In [ ]:
#describe the data
gdf.describe()
Out[ ]:
co2_std co2_sdev lat lon
count 13104.000000 13104.000000 13104.000000 13104.000000
mean 0.000002 0.000330 -0.005495 -1.250000
std 0.000001 0.000148 52.528319 103.924508
min 0.000000 0.000000 -90.000000 -180.000000
25% 0.000002 0.000394 -46.000000 -90.625000
50% 0.000003 0.000395 0.000000 -1.250000
75% 0.000003 0.000396 46.000000 88.125000
max 0.000005 0.000419 89.500000 177.500000
In [ ]:
gdf.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
Index: 13104 entries, 0 to 13103
Data columns (total 5 columns):
 #   Column    Non-Null Count  Dtype   
---  ------    --------------  -----   
 0   co2_std   13104 non-null  float64 
 1   co2_sdev  13104 non-null  float64 
 2   lat       13104 non-null  float64 
 3   lon       13104 non-null  float64 
 4   geometry  13104 non-null  geometry
dtypes: float64(4), geometry(1)
memory usage: 614.2 KB
In [ ]:
import plotly.express as px
fig = px.histogram(df_dim, x="co2_std", nbins=80, title="Histograma de co2_std")
fig.show()
In [ ]:
import plotly.express as px
fig = px.histogram(df_dim, x="co2_sdev", nbins=200, title="Histograma de co2_sdev")
fig.show()

Basado en el segso observado en los graficos, optamos por eliminar valores donde co2_std -> 0¶

In [ ]:
gdf2 = gdf.query("co2_std >= 0.000000001") #.plot(column='co2_std', cmap='OrRd', figsize=(15, 15), legend=True, scheme='quantiles')
In [ ]:
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.colors as mcolors
import matplotlib.cm as cm
In [ ]:
gdf2["co2_std"]
Out[ ]:
8        0.000003
91       0.000003
119      0.000004
144      0.000002
145      0.000002
           ...   
10939    0.000002
10940    0.000003
10941    0.000002
10942    0.000003
10943    0.000003
Name: co2_std, Length: 10803, dtype: float64

Distribucion de datos despues de eliminar valores iguales o cercanos a 0¶

In [ ]:
import plotly.express as px
fig = px.histogram(gdf2, x="co2_std", nbins=80, title="Histograma de co2_std")
fig.show()
In [ ]:
import plotly.express as px
fig = px.histogram(gdf2, x="co2_sdev", nbins=200, title="Histograma de co2_sdev")
fig.show()

Analisis de estacionariedad¶

In [ ]:
lati = gdf2["lat"]
loni = gdf2["lon"]
latlon = list(zip(lati, loni))
latlon = [str(x) for x in latlon]
fig,ax = plt.subplots()
plt.plot(latlon,gdf2["co2_std"],color='red')
plt.xticks(range(0,len(latlon),(round(len(latlon)/20))))
plt.xticks(rotation=60)
plt.xlabel("(Latitud, Longitud)")
plt.ylabel("Valor Co2_std")
plt.title("CO2 troposferico")
plt.show()
No description has been provided for this image
In [ ]:
lati = gdf2["lat"]
loni = gdf2["lon"]
latlon = list(zip(lati, loni))
latlon = [str(x) for x in latlon]
fig,ax = plt.subplots()
plt.plot(latlon,gdf2["co2_sdev"],color='purple')
plt.xticks(range(0,len(latlon),(round(len(latlon)/20))))
plt.xticks(rotation=60)
plt.xlabel("(Latitud, Longitud)")
plt.ylabel("Valor Co2_sdev")
plt.title("desviacion estandar CO2 troposferico")
plt.show()
No description has been provided for this image

Pruebas de estacionariedad¶

Para esto utilizaremos pruebas de raiz unitaria. Una raíz unitaria es una característica de los procesos que evolucionan a través del tiempo y que puede causar problemas en inferencia estadística en modelos de series de tiempo.

Un proceso estocástico lineal tiene una raíz unitaria si el valor de la raíz de la ecuación característica del proceso es igual a 1, por lo tanto tal proceso es no estacionario. Si las demás raíces de la ecuación característica se encuentran dentro del círculo unitario es decir, tienen un valor absoluto menor a uno. entonces la primera diferencia del proceso es estacionaria.

Algunas pruebas de raiz unitarias son:¶

Argumented dickey fuller¶

La prueba de Dickey-Fuller aumentada es una prueba de hipótesis. La hipótesis nula es que la serie de tiempo no es estacionaria, y la alternativa es que la serie es estacionaria. Por lo tanto, necesitamos encontrar un valor p lo suficientemente bajo como para rechazar nuestra hipótesis nula, lo que sugiere que la serie es estacionaria.

Phillips Perron¶

La prueba Phillips perron es una prueba de raiz unitaria basada en la prueba de dickey fuller y su misma hipotesis nula. y si el valor p es menor a la significanica rechazamos la hipotesis nula.

Kwiatkowski Phillips Schmidt Shin¶

Contrariamente a la mayoría de las pruebas de raíz unitaria , la presencia de una raíz unitaria no es la hipótesis nula sino la alternativa. Además, en la prueba KPSS, la ausencia de una raíz unitaria no es una prueba de estacionariedad sino, por diseño, de estacionariedad debil. Esta es una distinción importante ya que es posible que una serie de tiempo no sea estacionaria, no tenga raíz unitaria pero sea debilmente estacionaria. por lo tanto si el valor p es menor a la significancia se dice que la serie es debilmente estacionaria.

Dickey Fuller¶

La Prueba de Dickey-Fuller busca determinar la existencia o no de raíces unitarias en una serie de tiempo. La hipótesis nula de esta prueba es que existe una raíz unitaria en la serie. El Modelos de regresión puede ser escrito como:

$$ \nabla y_{t}=(\rho-1)y_{t-1}+u_{t}=\delta y_{t-1}+u_{t} $$

Donde ∇ es el operador de primera diferencia. Este modelo puede ser estimado y las pruebas para una raíz unitaria son equivalentes a pruebas δ = 0 (donde δ = ρ - 1). y al igual que en la prueba aumentada para rechazar la hipotesis nula el valor p obtenido de la prueba debe ser menor a nuestra significancia.

In [ ]:
# Librerias
import pandas as pd
import matplotlib.pyplot as plt 
from statsmodels.tsa.stattools import adfuller
from arch.unitroot import PhillipsPerron,KPSS,DFGLS,ADF
from prettytable import PrettyTable
In [ ]:
# ignore future warnings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
In [ ]:
gdf2.columns
Out[ ]:
Index(['co2_std', 'co2_sdev', 'lat', 'lon', 'geometry'], dtype='object')

Pruebas para CO2 Troposférico std¶

In [ ]:
serie = gdf2.iloc[:,0]
adf = ADF(serie)
pp = PhillipsPerron(serie)
kpss = KPSS(serie)
df = DFGLS(serie)
x = PrettyTable()
x.field_names = ["Unitary root test","Test Statistic", "P-value", "Critical Value for 5%", "Alternative Hypothesis:"]
x.add_row(["Argumented Dickey Fuller", adf.stat, adf.pvalue, adf.critical_values["5%"], adf.alternative_hypothesis])
x.add_row(["Phillips Perron", pp.stat, pp.pvalue, pp.critical_values["5%"], pp.alternative_hypothesis])
x.add_row(["Kwiatkowski-Phillips-Schmidt-Shin", kpss.stat, kpss.pvalue, kpss.critical_values["5%"], kpss.alternative_hypothesis])
x.add_row(["Dickey-Fuller GLS", df.stat, df.pvalue, df.critical_values["5%"], df.alternative_hypothesis])
x
Out[ ]:
Unitary root test Test Statistic P-value Critical Value for 5% Alternative Hypothesis:
Argumented Dickey Fuller -9.27058995597038 1.3301353820184625e-15 -2.8618085021011126 The process is weakly stationary.
Phillips Perron -101.07287313284054 0.0 -2.861807607138271 The process is weakly stationary.
Kwiatkowski-Phillips-Schmidt-Shin 3.5542876734971465 0.0001 0.4614 The process contains a unit root.
Dickey-Fuller GLS -8.486833870023165 6.730213808921115e-14 -1.9456491396867541 The process is weakly stationary.

En resumen, las pruebas de Dickey-Fuller Aumentada, Phillips-Perron y Dickey-Fuller GLS indican que el proceso es débilmente estacionario, ya que los valores P son muy cercanos a cero, y los estadísticos de prueba están por debajo de los valores críticos.

Por otro lado, la prueba KPSS indica que el proceso contiene una unidad raíz, lo que significa que no es estacionario. Esto crea una discrepancia entre las pruebas, y en este caso, se daría preferencia a las pruebas de Dickey-Fuller Aumentada, Phillips-Perron y Dickey-Fuller GLS, que indican estacionariedad. Sin embargo, es importante considerar la interpretación en el contexto de los datos y los objetivos del análisis. Esta incertidumbre debido a la discrepancia en los resultados, nos motivan a considerar otras técnicas de análisis.

Pruebas para la desviacion estandar del Co2¶

In [ ]:
serie = gdf2.iloc[:,1]
adf = ADF(serie)
pp = PhillipsPerron(serie)
kpss = KPSS(serie)
df = DFGLS(serie)
x = PrettyTable()
x.field_names = ["Unitary root test","Test Statistic", "P-value", "Critical Value for 5%", "Alternative Hypothesis:"]
x.add_row(["Argumented Dickey Fuller", adf.stat, adf.pvalue, adf.critical_values["5%"], adf.alternative_hypothesis])
x.add_row(["Phillips Perron", pp.stat, pp.pvalue, pp.critical_values["5%"], pp.alternative_hypothesis])
x.add_row(["Kwiatkowski-Phillips-Schmidt-Shin", kpss.stat, kpss.pvalue, kpss.critical_values["5%"], kpss.alternative_hypothesis])
x.add_row(["Dickey-Fuller GLS", df.stat, df.pvalue, df.critical_values["5%"], df.alternative_hypothesis])
x
Out[ ]:
Unitary root test Test Statistic P-value Critical Value for 5% Alternative Hypothesis:
Argumented Dickey Fuller -4.438444574142901 0.00025362816868106694 -2.8618085769516295 The process is weakly stationary.
Phillips Perron -25.853410981967464 0.0 -2.861807607138271 The process is weakly stationary.
Kwiatkowski-Phillips-Schmidt-Shin 12.415548071750424 0.0001 0.4614 The process contains a unit root.
Dickey-Fuller GLS -0.7686066800043696 0.3939327633407306 -1.945649700954751 The process is weakly stationary.

En este caso, la Prueba de Dickey-Fuller Aumentada y la Prueba de Phillips-Perron indican que el proceso es débilmente estacionario, ya que los valores P son muy bajos, y los estadísticos de prueba están por debajo de los valores críticos. Estas pruebas respaldan la estacionariedad de los datos.

Por otro lado, la Prueba de Kwiatkowski-Phillips-Schmidt-Shin (KPSS) sugiere que el proceso contiene una unidad raíz, lo que significa que no es estacionario. Esta prueba está en contradicción con las otras dos pruebas.

La Prueba de Dickey-Fuller GLS no es tan concluyente, ya que el valor P es relativamente alto, y el estadístico de prueba no rechaza la hipótesis de no estacionariedad, pero tampoco la confirma con firmeza.

En contexto¶

Dado el contexto, es importante considerar que la estacionariedad de estas series temporales puede ser crucial en análisis geoestadísticos relacionados con el CO2 atmosférico, ya que las tendencias y patrones temporales pueden influir en la interpretación de los datos y en la toma de decisiones.

Las pruebas de raíz unitaria que has realizado indican que la Prueba de Dickey-Fuller Aumentada y la Prueba de Phillips-Perron sugieren que las series de datos "co2_std" y "co2_sdev" son débilmente estacionarias, lo que implica que pueden ser utilizadas en análisis estadísticos con la suposición de estacionariedad.

Por otro lado, la Prueba de Kwiatkowski-Phillips-Schmidt-Shin (KPSS) sugiere que las series de datos contienen una unidad raíz, lo que implicaría no estacionariedad. Esta discrepancia entre las pruebas podría deberse a diferentes supuestos y enfoques de las pruebas.

En esta situación, consideramos las siguientes opciones:

Evaluar visualmente los datos: Realiza gráficos de las series temporales para identificar tendencias, estacionalidades u otros patrones visuales que puedan ayudarte a tomar una decisión más informada sobre la estacionariedad.

Realizar más pruebas: Podrías considerar realizar pruebas adicionales o explorar técnicas estadísticas más avanzadas para evaluar la estacionariedad de las series de datos.

Diferenciación: Si no puedes llegar a una conclusión clara sobre la estacionariedad de las series de datos, podrías aplicar técnicas de diferenciación para transformar los datos y hacerlos estacionarios antes de continuar con tu análisis geoestadístico.

Logaritmos: En nuestro caso particular y debido que los valores son pequeños, aplicar logaritmos aún puede ser una estrategia útil para estabilizar la varianza y reducir la magnitud de las fluctuaciones en la serie temporal. Sin embargo, es importante tener en cuenta que, en este caso, el logaritmo natural (ln) puede ser más apropiado, ya que tiende a funcionar mejor con valores pequeños. Como datos son muy pequeños, es probable que tengan una magnitud similar a la de los valores logaritmicos resultantes. Esto puede facilitar la interpretación de los resultados y hacer que los datos transformados sean más adecuados para análisis groestadísticos, incluyendo pruebas de estacionariedad.

In [ ]:
kpss.alternative_hypothesis
Out[ ]:
'The process contains a unit root.'
In [ ]:
kpss.summary()
Out[ ]:
KPSS Stationarity Test Results
Test Statistic 3.554
P-value 0.000
Lags 58


Trend: Constant
Critical Values: 0.74 (1%), 0.46 (5%), 0.35 (10%)
Null Hypothesis: The process is weakly stationary.
Alternative Hypothesis: The process contains a unit root.
In [ ]:
# Crea una figura de Matplotlib
fig, ax = plt.subplots(figsize=(12, 8))

# Dibuja el mapa del mundo de GeoPandas como fondo
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.boundary.plot(ax=ax, linewidth=1, color='black')

# Utiliza Seaborn para dibujar los puntos y personalizar la paleta de colores
sns.scatterplot(ax=ax,x= 'lon', y= 'lat', data = gdf2, hue = 'co2_std', palette='coolwarm')

# Agrega opacidad a los puntos para diferenciar la concentración
alpha = 0.3
ax.collections[0].set_alpha(alpha)

# Personaliza el aspecto del mapa
ax.set_title('Distribución de CO2 en la Troposfera Libre con Opacidad')
plt.axis('off')  # Desactiva los ejes

sm = cm.ScalarMappable(cmap='coolwarm')
sm.set_norm(mcolors.Normalize(vmin=min(gdf2["co2_std"]), vmax=max(gdf2["co2_std"])))
# Agrega una barra de colores (colorbar)
cbar = plt.colorbar(sm, ax=ax)
cbar.set_label('Concentración de CO2 (co2_std)')


# Muestra el mapa con los datos
plt.show()

/tmp/ipykernel_2109/3664182266.py:5: FutureWarning:

The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
No description has been provided for this image
In [ ]:
# Crea una figura de Matplotlib
fig, ax = plt.subplots(figsize=(12, 8))

# Dibuja el mapa del mundo de GeoPandas como fondo
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.boundary.plot(ax=ax, linewidth=1, color='black')

# Utiliza Seaborn para dibujar los puntos y personalizar la paleta de colores
sns.scatterplot(ax=ax,x= 'lon', y= 'lat', data = gdf2, hue = 'co2_sdev', palette='viridis')

# Agrega opacidad a los puntos para diferenciar la concentración
alpha = 0.3
ax.collections[0].set_alpha(alpha)

# Personaliza el aspecto del mapa
ax.set_title('Distribución de CO2_sdev en la Troposfera Libre con Opacidad')
plt.axis('off')  # Desactiva los ejes

sm = cm.ScalarMappable(cmap='viridis')
sm.set_norm(mcolors.Normalize(vmin=min(gdf2["co2_sdev"]), vmax=max(gdf2["co2_sdev"])))
# Agrega una barra de colores (colorbar)
cbar = plt.colorbar(sm, ax=ax)
cbar.set_label('Concentración de CO2 (co2_sdev)')


# Muestra el mapa con los datos
plt.show()

/tmp/ipykernel_2109/1589977332.py:5: FutureWarning:

The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
No description has been provided for this image
In [ ]:
gdf2.columns
Out[ ]:
Index(['co2_std', 'co2_sdev', 'lat', 'lon', 'geometry'], dtype='object')
In [ ]:
import numpy as np
import pandas as pd
from geostatsmodels import utilities, variograms, model, kriging, geoplot

# Extrae las indices y los valores de "co2_std" del GeoDataFrame
indices = gdf2.geometry.apply(lambda geom: [geom.x, geom.y]).to_list()
values = gdf2[["lon","lat" , "co2_std"]]
In [ ]:
P = values.to_numpy()
P
Out[ ]:
array([[-1.60000000e+02,  8.95000000e+01,  3.03339266e-06],
       [ 4.75000000e+01,  8.95000000e+01,  3.42648696e-06],
       [ 1.17500000e+02,  8.95000000e+01,  3.83136512e-06],
       ...,
       [ 1.72500000e+02, -6.00000000e+01,  2.16231664e-06],
       [ 1.75000000e+02, -6.00000000e+01,  2.50185760e-06],
       [ 1.77500000e+02, -6.00000000e+01,  2.77377723e-06]])
In [ ]:
pw = utilities.pairwise(P)
geoplot.hscattergram(P,pw,10,5)
geoplot.hscattergram(P,pw,20,5)
geoplot.hscattergram(P,pw,30,5)
In [ ]:
gdf["co2_std"].mean()
In [ ]:
# Calcula el semivariograma

tolerance = gdf["co2_std"].std()
lags = np.arange( tolerance, gdf["co2_std"].mean(), tolerance * 2 )
sill = np.var(P[:,2])
svm = model.semivariance( model.spherical, [ 4000, sill ] )
geoplot.semivariogram( P, lags, tolerance, model=svm )

You have more than 10,000 data points, this might take a minute.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/workspaces/Geo/note.ipynb Celda 25 line 7
      <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X30sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a> sill = np.var(P[:,2])
      <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X30sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a> svm = model.semivariance( model.spherical, [ 4000, sill ] )
----> <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X30sdnNjb2RlLXJlbW90ZQ%3D%3D?line=6'>7</a> geoplot.semivariogram( P, lags, tolerance, model=svm )

File ~/.python/current/lib/python3.10/site-packages/geostatsmodels/geoplot.py:87, in semivariogram(data, lags, tol, model)
     75 '''
     76 Input:  (data)    NumPy array with three columns, the first two
     77                   columns should be the x and y coordinates, and
   (...)
     84 Output:           empirical semivariogram
     85 '''
     86 # h, sv = variograms.semivariogram(data, lags, tol)
---> 87 vdata = variograms.semivariogram(data, lags, tol)
     88 h, sv = vdata[0], vdata[1]
     89 sill = np.var(data[:, 2])

File ~/.python/current/lib/python3.10/site-packages/geostatsmodels/variograms.py:71, in semivariogram(data, lags, tol)
     63 def semivariogram(data, lags, tol):
     64     '''
     65     Input:  (data) NumPy array where the first two columns
     66                    are the spatial coordinates, x and y
   (...)
     69     Output: (sv)   <2xN> NumPy array of lags and semivariogram values
     70     '''
---> 71     return variogram(data, lags, tol, 'semivariogram')

File ~/.python/current/lib/python3.10/site-packages/geostatsmodels/variograms.py:121, in variogram(data, lags, tol, method)
    119     v = [covariance(data, indices) for indices in index]
    120 # bundle the semivariogram values with their lags
--> 121 return np.c_[lags, v].T

File ~/.local/lib/python3.10/site-packages/numpy/lib/index_tricks.py:418, in AxisConcatenator.__getitem__(self, key)
    414     # concatenate could do cast, but that can be overriden:
    415     objs = [array(obj, copy=False, subok=True,
    416                   ndmin=ndmin, dtype=final_dtype) for obj in objs]
--> 418 res = self.concatenate(tuple(objs), axis=axis)
    420 if matrix:
    421     oldndim = res.ndim

ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 0, the array at index 0 has size 1 and the array at index 1 has size 0
In [ ]:
# Ajusta un modelo de semivariograma
# Puedes elegir diferentes modelos (spherical, exponential, etc.) y ajustar sus parámetros.
vgm = model.vgm(gdf2, "co2_std", "exponential", 10, 10)

# Grafica el semivariograma y el modelo ajustado
import matplotlib.pyplot as plt

plt.scatter(svm['lag'], svm['val'], label='Semivariogram', s=60)
plt.plot(vgm['distance'], vgm['gamma'], label='Model', color='red')
plt.xlabel('Lag Distance')
plt.ylabel('Semivariance')
plt.legend()
plt.show()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/workspaces/Geo/note.ipynb Celda 26 line 3
      <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a> # Ajusta un modelo de semivariograma
      <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=1'>2</a> # Puedes elegir diferentes modelos (spherical, exponential, etc.) y ajustar sus parámetros.
----> <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=2'>3</a> vgm = model.vgm(gdf2, "co2_std", "exponential", 10, 10)
      <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a> # Grafica el semivariograma y el modelo ajustado
      <a href='vscode-notebook-cell://codespaces%2Bprobable-trout-464xww99w6vhj67r/workspaces/Geo/note.ipynb#X34sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a> import matplotlib.pyplot as plt

NameError: name 'model' is not defined